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I. 


INTRODUCTION 


In recent years, significant activities have been undertaken 
to improve the efficiency, stability and range of application of 
Newton Raphson (NR) type solution schemes in handling nonlinear 
finite element (FE) simulations. This has led to the development 
of hyper-linear , circular elliptic and piecewise 

continuous ^ constraint surfaces. These are used to control the 
I NR iterative process, thereby, enabling the handling of transition 
in Jacobian definiteness, as well as, a wide variety of geometric 
and material nonlinearities Such work has emphasized the 

following developmental aspects namely: 

i) Streamline Jacobian (stiffness matrix) updating for 
instance as in the use of the BFGS scheme ^8-10], 

ii) Introduce use of constraints to bound iterations associated 

with the incremental Newton Raphson (INR) algorithm, thus, 
enabling the handling of problem with indefinite Jacobian 
properties and 

iii) Introduce self adaptive attri butes enabl ing automatic load 
stepping 

The above work, in particular the constrained methodology has 
significantly extended the range of INR based schemes, regardless several diffi- 
culties are still present. In addition to the storage intensiveness of 
the scheme, the user is faced with the lack of apriori or automatic 
procedures enabling adjustment of size, shape and orientation of the 
constraint surface, so as to enhance numerical robustness. In this con- 
text, the user must undergo the trial and error process of defining 
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load stepping. Depending on the constraint function employed, this 
involves the specification of the various geometric attributes of the 
bounding surface. For instance, in the case of the hyperelliptic con- 
straint surface (HECS), the aspect ratio and overall ordinate and 
abscissa bounds must be established on a global level. 

To circumvent the foregoing difficulties, this paper develops 
numerical strategies which enable the automatic adjustment of the 
size, shape and orientation of the constraint surface associated with 
full (FINR), modified (MINR) and BFGS updated Newton Raphson FE equation 
solvers. This will be achieved through the introduction of a hierarchy 
of localized constraints which enable the control of the iterative ex- 
cursions of various levels of the governing field variables. The approach 
taken is general in that a wide variety of constraint surfaces can be 
handled by the methodology. In addition to the general development, the 
results of several extensive benchmark activities is also presented. Here, 
special emphasis is given to illustrating the flexibility and robustness 
of the locally bound approach to controlling successive globally con- 
strained INR type solution algorithms. In the sections which follow, 
detailed discussions will be given on the nature and extent of shortcomings 
of constrained approaches, the development of reshaping and resizing 
schemes from a 1-D point of view, the development of local constraint 
methodologies via FE generalization. Finally, the results of extensive 
benchmark examples demonstrating enhanced numerical characteristics will 
be presented. 

1 1 , CONSTRAINED NEWTON RAPHSON SCHEMES: STRENGTHS/SHORTCOMINGS 


As noted earlier, INR type solution schemes involving either 

modified, straight or BFGS type stiffness updating suffer from several 
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very basic shortcomings To simplify the discussion and develop- 

ment, the INR is reinterpreted to be a linearly bound hyperl inear ex- 
trapolation of the solution curve via a full, modified or BFGS updated 
stiffness. Specifically, noting Fig. 2-1, the succession of extrapola- 
tion are bound by the intersection with the hyperplane defined by the 
designated load level. For such a process, it is difficult to guarantee 
the requisite intersections. This is especially true for structural 
stability problems involving buckling with turning point behavior where 
slope indefiniteness may occur. 

To bypass such problems, the surface bounding the INR extrapolation 
can be reinterpreted to be a hyperplane which is oblique to F = constant 
or alternatively a different function. Such an interpretation naturally 
yields the approaches taken by Wempner Riks Crisfield ^ and 
Padovan, et al. Specifically, noting Fig. 2-2, this has led to 

the use of constant arc lengths oblique hyperlines, and circular 

elliptic hyperbolic as well as, piecewise continuous ^ 

bounding constraint surfaces. While all such methodologies can handle 
the turning point problems associated with structural stability, as shown 
in Fig. 2 -3, those which involve closed constraint curves more readily 
guarantee the requisite intersection with the INR extrapolation 

In spite of the advanced features of the foregoing schemes, the con- 
strained approach also suffers from several very important fundamental 
shortcomings. While various of these are still outgrowths of the metho- 
dologies INR base, several are a direct function of the use of the con- 
straint concept. Overall, the shortcomings are defined by the following: 

i) Unlike the F = constant constraint surface, all the other 

surfaces yield intersections with the INR extrapolations which , 
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due to system nonlinearity remain apriori unknown; this makes 
the calculation for specific so called target loads awkward. 

ii) While the various constraint surfaces self adaptively adjust 
load stepping during successive iterations, such incrementation 
is a direct outgrowth of their intrinsic geometries. Hence, 
the user is still left with the oftentimes puzzling problem of 
defining their requisite shape, aspect ratio, obi i queness ,*arc 
length , radi i , etc . 

iii) Because of the lack of apriori knowledge of the intersection of 
the constraint and I NR extrapolation, it is difficult to handle 
mixed problems which may involve some combination of prescribed 
displacement and loading requirements or targets. 

iv) Oversized constraint surfaces may initiate a nonintersection in 
hyperspace and hence solution failure. (See Fig. 2-4). In con- 
trast, an overly tight surface can result in an inefficient and 
hence, expensive iteration process. 

v) Generally, since the intersections of the constraint curve and 

solution surface are not apriori known, typically the constrained 
methodology will converge to points significantly reduced from 
expected load steps: Noting Fig. 2-5, this leads to under uti- 

lized load steps and consequently numerical inefficiency. 

* 

In the following sections, several methodologies will be imple- 
mented to bypass the foregoing shortcomings. As will be seen, these 
are general enough to be employed with most of the constraints in 
use today. 
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HI. IMPROVED CONSTRAINT SCHEME: 1-D DEVELOPMENT 


To simplify the development of the improved constraint methodology, 
we shall first consider modifications to the solution of 1-D nonlinear 
equations. Recalling the comments of the preceding section, it follows 
that the shortcomings of the NR based constraint schemes fall into 
several main categories namely: 

i) The need to develop methodologies which - 

• Enable convergence to specific load and/or 
deflection states and; 

• Allow full use of a given load step; 

ii) Develop schemes which automate - 

• Sizing of constraint surface and; 

• Reshape and reorient constraint surface 

so as to improve convergence characteristics. 

To circumvent the foregoing shortcomings, several different 
modifications will be introduced into the constrained methodology. 

The main thrust of these will be twofold, namely: 

i) To provide for self adaptive modifications of constraint 
geometry enabling convergence to designated target load 
and/or deflection state and; 

ii) To employ localized constraints so as to enable essentially 
automatic sizing of constraint while still maintaining 
numerically stable convergence characteristics. 

While the geometries of the various candidate constraint hyper- 
surfaces are unique, there are basic coordinate transformations which 

will enable flexible reshaping and reorientation. Specifically, for 
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this paper, we shall consider the use of translations, rotations 
and coordinate stretches. Obviously, dependent on the type of con- 
straint geometry used, some of these schemes are more advantageous 
than others. For example, consider the use of the oblique hyperplane, 

Fig. 3-1 illustrates combinations of rotation and rotation-translations. 
While both schemes possess the requisite targeting capabilities, the 
piecewise continuous rotation-translation scheme may yield a more 
stable numerical convergence process. Note, in the vicinity of turning 
points, the rotation process of both methodologies would have to be 
stopped to guarantee an intersection with the solution curve. 

Considering the piecewise hyperbolic constraint. Fig. 3-2 de- 
monstrates the use of translations to enable requisite targeting. As 
with the oblique hyperplanar constraint, in the vicinity of turning 
points, the translation process is stooped to insure the requisite 
intersections . 

For the elliptic constraint, abscissa or ordinate coordinate 
stretches can be used to enable the requisite targeting. Such a pro- 
cess is illustrated in Fig. 3-3. In particular, the iterations depicted 
involve an abscissa stretched hyperelliptic constraint surface AHECS. 

To establish the requisite algorithmic formulation, we shall consider 
the following 1-0 nonlinear equation namely: 

G(Y) = F (3.1) 

The NR scheme associated with (3.1) is given by the expression 

Vl ‘ Y n + < 5 - G < Y n» /G '< V n> i G '<V = 37 (G(Y n>> < 3 - 2 > 

such that $ is the designated value of F, In terms of the constraint 
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concept, the iterations defined by (3.2) can be bound by rescaling $ 
via the parameter x. The choice of X is given by the intersection of 
the NR extrapolation and the constraint curve. For example, such inter- 
sections are defined in Fig. 2-2. Overall, such a process is given by 
the following relations: 

i) Constrained NR algorithm; 

Vl ■ Y n + <Vl 4 - G(Y n »/S'(Y„) (S.3) 

i i ) Constraint curve ; 

C (Y , *) = 0 (3.4) 

Such that the intersection takes the form 


I (x 


n+1 


Y n+1 ^ 


0 . 


To determine x^ + ^ , the specific form of C must be chosen. In 
the case of oblique linear, circular, elliptic and piecewise hyperbolic 
curves/surfaces, we have that 

i ) Obi ique 1 inear 

F = a Y + » (3.5) 

ii) Circular 

F 2 + Y 2 = $ 2 (3.6) 

ili) Elliptic 

F 2 + yY 2 = * 2 (3.7) 

iv) Piecewise Hyperbolic 

(F - *)(Y - Y ) = u (3.8) 
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where the parameters 3, *>, y and p are user defined. As noted earlier, 
such relationships can be modified to admit translation rotation and 
coordinate stretching. 

For instance, considering the oblique linear constraint, rotation 
can be introduced by varying the selection of 3 in the range 3 e [ 0 , u); 
p < ® for F > 0 and 3 e [--j, 0); u < » for F < 0. The actual variation 
of 3 can be undertaken in several different ways, namely 

i) Angular increments, Fig. 3-4; 

ii) Force increments, Fig. 3-5; 

iii) MNR estimated aY, Fig. 3-6; 

iv) BFGS-NR estimated aY, and via 

v) BFGS-NR or MNR estimated aF. 

Note in the case of fixed aY incrementation, the sweeping of the 
oblique linear constraint is modified to handle iterative steps in the 
vicinity of the designated load and/or deflection state. Specifically, 
the constraint is set to F = t. Overall, the algorithm associated with 
the fixed MNR estimated aY sweep of the oblique linear constraint requires 
a knowledge of successive intersection points. That is, considering the 
(n+l) th rotation, the associated MNR extrapolation and oblique linear 
constraint curve are given by the following relations: 

Y = Y n + (F - G(Y n ))/G' (3.9) 


and 


F = 3 n+1 Y ♦ * 


(3.10) 
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In terms of (3.9) and (3.10), the intersection I defined in Fiq. 3-7 

n 3 

is given by 


Y n = Y n + <Vl 0 - G <V> /G ' 


(3.11) 


Vi 5 = * 


., Y + 
n+1 n 


(3.12) 


Since 


is given 

and Vi 


* iV (3.13) 

(3.11) and (3.12) can be used to define the unknowns a .j 
After several manipulations we have that 


a n+l 


ts - " — “ (G 1 AY + G(YJ - $} 

Y + iY n 

n 


(3.14) 


A* , = 5-^ — (G 1 ay + G(Y°) - $} (3.15) 

n 1 $(Y + &Y) n 

n 

Note, for the BFGS or fully updated schemes, G' appearing in (3.14) and 
(3.15) can be replaced by its corrected values namely G 1 (Y n ). 

The swept oblique linearly constrained MNR algorithm has two basic 
forms. Namely for Vi « ’■ (3.14) and (3.15) are used to perform the 

l 

necessary iterations. When X n+ -j > 1, then sweeping is stopped and iteration 
continues as per (3.2). Note, in the vicinity of turning or buckling points, 
the foregoing algorithm must be modified to insure the requisite intersection. 
This can be achieved in several ways that is: 

i) The curvature can be monitored as per the G* (Y ) parameter; 
when below critical level sweeping is stopped; 
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ii) Alternatively, a monotonicity check of the G(Y n ) sequence can 
be used to terminate sweeping. 

In a similar context, the alternative constraint curves can be 
modified to include rotations, translations and stretches. For example, 
we .shall consider the more versatile elliptic surface defined in Figs. 
2-5 and 3-3, Considering the (n+l)^ h successive stretch depicted in 

I 

Fia. 3-8, the following relations define the intersection I namely 

n J 

<Vi 5 - G(V > G ’ (3 - 16) 


<*; +1 »> 2 + ( vj 2 = « 2 


n+l v n‘ 


(3.17) 


»’ « 

Due to (3.13), solving (3.16) and (3.17) for X n+ ^ and u n+ -| we have that 


Vi ■? !G ': ^ + G(Y n> > 


(3.18) 


u n+ 1 = — 5 —^ — ■=-* {$ 2 - (G AY - G(Y )) } 

n 1 (Y + AY) Z n 


(3.19) 


Like the oblique linear constraint, the elliptically bound NR 

can be employed in a variety of ways. These include: 

i) After each stretch of the elliptic constraint curve, several 
iterations can be used to stabilize the numerical convergence 
process; 

ii) The ongoing stretch of the elliptic constraint can be used to 
delimit successive INR type iterations; 

iii) The constraint can be employed with either BFGS, full or 
modified slope updating; 
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z 

iv) For x p - 0 (1), lY can be set to » (large number) to enable 
convergence to designated load level; 

l 

v) When the definiteness of G varies as in the vicinity of turning 
points, the stretching process is terminated. 

As noted by Padovan and Arechaga so-called safety zones can 
be established for the constrained scheme wherein convergence can be 
guaranteed under the appropriate conditions. For instance, let us con- 
sider the case where successive iterations are undertaken between abscissa 
or ordinate stretches of the ellipse. Recalling (3.16) and (3.17), once 
- n +l is fixed by the stre tch, (3.19), we must insure that subsequent 
iterations lie in the appropriate safety zone. This is determined by 
monitoring the intersection points of subsequent iterations. Specifically, 
noting Fig. 3-8, it follows that the £ tb intersection of the n tb ordinate 
stretch is given by 

V n * C’ + < x n+l 5 ' SOj'hVG' (3-20) 


( x n+l ” 2 + “n+l 


(3.21) 


Solving for x n+ , we obtain the following quadratic identity, that is 

£ 


<Cv> 2 C Un + 2 q+1 C 2nn + C 3tn * 0 


where 


C Un = * 2 * < i/G '> 2 “n+l 


(3.22) 


(3.23) 


'2in 


" (C 1 ' 6(y‘- 1 )/G , )(*/g‘) u n+ , 


(3.24) 
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(3.25) 


C = (Y £_1 
L 3tn u n 


Giv'-'j/s ') 2 vi 


- 4> 


In terms of (3.22), we yield the expression 


' n+1 


1 

7 


C Un 


^ C 2m^ " C Un C 3tn^ } 


(3.26) 


Following the formalism of Padovan and Arechaga ^ , for the 

current purposes, the safety zone is defined by those values of Y and 
F where the discriminant of (3,26) is positive definite. Noting (3.23) - 
(3.24), such a requirement yields that 


^ C 2tn^ - C Un C 3tn 

(Y* +1 - G(Yj-V) 2 (*/<*' ) 2 ^ n +l )2 - 


U 2 - (o/G 1 ) 2 Vl> '(C 1 - )/g’) z „ n+1 - * 2 ) (3.27) 

Such a criterion enables the appropriate choice of &Y and hence the 
sizing of the ellipse. 

The foregoing discussion has developed a variety of methodologies 
enabling the convergence to specified force and/or deflection states. 
Generally, this has involved the reshaping or sizing of the constraint 
surface. In the section to follow, the scheme will be automated through 
the use of localized constraints which retain targeting capabilities, 
while also enabling generalization to large scale multidegree of free- 
dom simulations. 
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IV. EXTENSIONS TO N-D: FE DEVELOPMENT 


To automate the control of the resizing/reshaping of the con- 
straint surface for multidimensional problems, the use of localized 
bounds will be introduced. This will be facilitated by employing a 
nonlinear FE formulation. Specifically, assuming that the standard 
2 nc * Piola Kirchhoff stress S^ and Lagrangian strain ^3 
tensor combination are employed, the governing field equations take 
the form 


Hj ( V*ik*^» + f o) -o 


S i j " S ij ^hr L 22 ’ 0 • 


(4.1) 

(4.2) 


such that 


[13] 


s . . = !sl !!i !!i 0 

P 3X 3X Sot 

01 5 


(4.3) 



(u 


i.J 


J.i 




(4.4) 


with p Q , p defining the initial and current densities, Euler 
coordinates, a^ Lagrangian coordinates, a., the Cauchy stress and 
u. are the components of the displacement vector. For nonlinear 
materials, generally (4.2) is cast in an incremental form. In this 
context 


iS ij ■ T ijik 1L tk 


(4.5) 


such that y. is the tangent stiffness of the medium. 

Using the usual displacement type formulation, the displacement 
vector of a given element can be related to its associated nodal de- 
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, fi4l 

formations via the expression L J 

u = [N(x)] Y (4.6) 

such that 

u * (u-j » U 3 ) (4.7) 

and [N], Y respectively denote the shape function and nodal deflection 

column vector : Employing (4.1 - 4.7) inconjunction with the virtual 

n 41 

work principle, the following FE formulations are obtained namely L J 

/ c (L) S dv 
R 


f [B*] T S dv * F (4.8) 

where Fisa column vector defining nodal forces, 6( ) the variational 
operator and 

(S) T = (S-| 1 , $ 22 ’ S 33> s i 2 * S 23’ S 13^ ^ 4 ' 9 ^ 

«(L) * [B*] «(Y) (4.10) 

with ^ 

[B*] = [B] + [B n (Y)] [G] (4.11) 

(L) - (i -11 > ^ 22 ' ^*33’ ^12* *"23 ’ ^"13^ (^*i2) 


The straight INR algorithm typically used to solve (4.8) is given 
by the expression El4,15] 


(jJiW ' Jm+l 


(4.13) 


where k defines the load step, m the iteration count for the given 
step and 
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(4.14) 


-m+l “ § (fj 

F k = F k_ ’ + iF k 

[K(V k )] ■ S {[G] T [S(Y k )][G] 

R 


+ tM_Y*))V(Y k )][B*(Y k )] )d, (4.16) 

k u 

such that [S(Y )] is the prestress matrix defined at Y k and [D(Y k )l 
~ UI ~m L J 

is the tangent material stiffness matrix namely 



■ CdCy5»3Cb* (y^)3 4Y k 


(4.17) 


For the MINR, the tangent stiffness [K ( Y)] is updated intermittently. 

For the constrained INR, the AY k + ^ excursion is bound by global 
level surfaces cast in terms of normed field variables. Specifically, 
considering the HECS, we have from Fig. 4-1 that 

II fll 2 + v II y || 2 = || A F k || 2 (4.18) 


where the Euclidean norm takes the form 

II y II 2 = E (y i ) 2 (4.19) 

i 

Typically, the choice of y and the target load is user defined. As 
can be seen from (4.18), such variables represent global level parameters 
and hence, are difficult to select. In as much as normed constraint 
equations like (4.18) or its linear; circular, or piecewise hyperbolic 
counterparts represent mathematically convenient expressions, they do not 
provide for control at local levels. This follows from the fact that 
norms represent gross upper bounds on the various excursions. In this 
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context they provide neither information on nor control of individual 
degrees of freedom. This situation applies to all the various types 
of field variables namely the deflections, deformation gradients, 
stresses, strains, strain energy density, temperature, thermal gra- 
dients velocities, accelerations etc. 

The foregoing difficulties can be partially circumvented by par- 
titioning the governing dependent finite element vectorial field.vari- 
ables into a hierarchy of different groups. At the lowest position 
of the hierarchy are individual field variables. These can be located 
either at given nodes or integration points. The overall levels which 
such groups can be partitioned consists of: 

i) Individual ^degrees of freedom; 

ii) Vectorial and tensorial quantities at nodes and 
integration points; 

iii) Whole elements and substructufal regions; 

iv) Locally averaged quantities as in nodal or integration 
point extrapolation-interpolations; and 

v) Material and non! inear groupings. 

Based on such a partitioning concept, the global norm of say Y is 
itself a summation of the various hierarchy of localized norms. In 
particular, || Y]| can be rearranged as follows: 

(II Y|| ) 2 = Z (Y.) 2 
i 

= E (il Y,j[ ) 2 (4.20) 

l 

where the hierarchy of vectors Y fc ;t =1,2,. . . define the various 
groups within the global vector Y namely 
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y t = {y|, f 2 , yJ, . . .} 

Similar hierarchies of groups can be established for all the de- 
pendent field variables. 

Overall, (4.20) can be applied at two levels namely, globally 
and locally. At the global level, ! j Y I j can be both employed within a 
constraint used to say control the overall iteration process of the 
INR so as to enable handling of Jacobian indefiniteness, as well as, 
within convergence tests. On the local scale, the various groups Y 
and their associated norms (| Y |] can be employed to delimit excursions 
generated during the overall iteration process. Such a dual level 
approach enables a higher degree of solution control. As will be 
seen later, this enables the use of much larger load/deflection incre- 
ments while still maintaining good convergence characteristics. 

The foregoing partitioned norms can be used to establish localized 
constraints which scale the load stepping. Like the global level app- 
lication, the local constraint can be chosen from a variety of function 
types namely constant, linear, circular, elliptic, piecewise hyperbolic 
or perhaps more general bounding functions. The manner of application 
depends somewhat on the field variable excursion being controlled. Re- 
gardless, several types of applications can be outlined. Overall, these 
can be categorized into two main types namely 

i) Greatest local upper bounds; and 

ii) Intermediate local bounds. 

For both procedures, excursions in the various field variables are first 

f8i 

obtained via the standard INR scheme with either full modified or BFGS L J 
type updating. Note, in the case of greatest upper bound schemes, the 
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constraint is chosen so that successive iterations converge essen- 
tially to the bound set. For the intermediate scheme, the constraint 
is chosen so that convergence is directed to limit points in the 
interval set by the bound. 

To illustrate such approaches we shall first consider the use of 
displacement type controls. In terms of successive iterate excursions, 
a local upper bound test (LUBT) is required to gage the various par- 
titioned groups. Employing the usual Euclidian norm to the partitioned 
groups yields the following LUBT, that is 


!1 AYJ| < u[l , u] (4.21 ) 

such that a^ defines the group upper bound. Its choice will be 
discussed later. If (4.21) is satisfied, then no rescaling of the 
load increment is necessary. In cases where (4.21) is violated, then 
a one parameter stret.ch/contraction is applied to the load step. In 
particular, (4.14) is recast in the form 

Ci = Ci C' 1 + C - ' s <?L )dv> < 4 - 22 > 


where 


min (- 


aY 


); 


II <0 II 

11 mi," £ 


ms , 1 


m+1 


= 




n II AY^||< A 


(4.23) 


It [l.u] 


Note the upper bound constraint defined by (4.22 and 4.23) can 
either be continuously or intermittently updated. If intermittent 
updating occurs, then the family is fixed at the beginning of each 


load increment. 
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For intermediate type bounds, condition (4.21) can be employed to 
resize the global constraint surface. Such a process generally leads 
to a succession of iterate excursions which are "smaller" or more con- 
servative than those obtained by the previous upper bound approach. 

As an example, consider the scaling of the abscissa dimension of say 
the HECS. This is achieved by adjusting the parameter y which controls 
the aspect ratio. Specifically, based on (4.12), the initial iterative 
excursion associated with say the k^ increment takes the form 



• [K(_Y*)) », k 


(4.24) 


Now requiring that (4.21) be satisfied for all partitions associated 
with Y, we obtain the following one parameter scaling of the HECS 
abscissa namely 


AF 


k ii 2 


e 


1 lUtfll 2 


(4.25) 


where 


r 


min ( ?. ■ iT".:— ); II Y ull > A 


Y k 

lu 




v. IS II A Yj t | 


> A. 


(4.26) 


(1 , y] 


Note for the modified abscissa rescaled HECS namely the AHECS, the 
constrained surface is updated at the beginning of each load step and 
is thereafter fixed during successive iteration steps associated with 
the given load increment. 

To enable the appropriate load targeting, y must be updated 
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during the iterative process. In this case (4.25 and 4.26) take the 
form 


k _ 
y m+l 


!| 4F 


"I 1 iv k /8 k 

r -r r 


2 


(4.27) 


such that 


11 L Im+U !| 


k _/ 

Vn 


vU II A Vi i !!-<** 


(4.28) 


Based on the foregoing scaling of y, the abscissa dinension of the 
AHECS is independent of the load step. In this context, under local 
constraint control, successive iterative excursions of the deflection 
field are largely independent of load step size so long as very large 
aF increments are employed. Specifically, noting Fig. 4.2, as load 
increment size is increased, the hyper tangent plane associated with 
the AHECS is essentially vertical. In this context, the intersections 
of the INR extrapolation and the AHECS occur near the upper bound con- 
straint on the displacement excursion. As will be seen from the forth- 
coming benchmark examples, the foregoing configurational properties of 
the AHECS enable a high degree of displacement targeting control, as well 
as, make the scheme largely load step size insensitive. 

In the context of such properties, the AHECS with localized dis- 
placement constraint enables the solution to progress to the requisite 
loading target under controlled upper bounded iterative increments. Such 
a load stepping process is more natural since the anticipated deflection 
in the structure can easily be evaluated directly from design level blue- 
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print considerations which define clearances, fits, acceptable struc- 
tural deformations etc. 


Pursuing this line of thought further, since strain is functionally 

9u. 

dependent on the various displacement gradients — l , it also seems 

aaj 

natural to employ localized constraints on admissible slope excursions. 

rU. 

In terms of the FE formulation, it follows that the various — - can be 

3a j 


recast in matrix format, that is 




777 ([N]V); j = 1, 2, 3 


(4.29) 


Based on (4.29), placing constraints on admissible 
following pointwise bound on Y, that is 



leads to the 



(Mr 1 




77 (M)' 1 JL (f) II 
J J 


< 


T*- (M ) 
° x j 


B <!> II 


5 

3X . 


([No- 


where r. defines the upper bound on 

J 

point of view, we have that 



(4.30) 

From an incremental 


II AY || < \\-£- ([N])' 1 || IT. (4.31) 

3 

for V x e R. If (4.31) is averaged say over a given element, then we 
obtain the more usable expression 
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(4.32) 


(‘Y k ) e ' 
Km+l' 1 


3X j 


-1 


dv 


i Li 


j * 1 , 2 , 3 

where V g is the volume of the e th element and (aY k +1 ) e the associated 
incremental nodal deflection. As can be seen, (4.32) represents a 
element level partitioning of the bounding process. 

In the context of (4.32), similar bounds can be 'defined for'strains, 
stresses as well as strain energy excursions. For instance, recalling 
the definition of strain, in terms of (4.4, 4.6, 4.11) we, see that 

■ (w ♦[»„(»$,, ms]) *i +1 < 4 - 33 > 

Considering say the e th element, the norm of (4.33) leads to the following 

l< 

bounds for admissible AY m+ -| namely 

^♦1 ■ <[B] + -L^ tl , < 4 -34) 


ii 'in+i ii * ii ««♦£»„<&, orar 1 <iii 

<ji([B] +[B n (^ +1 )][G])- 1 !| ||^ t1 !| . 

<11 ([B] + CB n (^ +1 )][G]>' 1 II B <4.35) 

k 

where 5 defines the upper bound on j| £L m+ -j || . Averaging over say the 
e*^ element yields the more useable constraint 

II < 4 Im+l I 6 " - V ell ' [B * ] dv H H e (4 ' 36> 

e 

Noting the relation for stress increments, the following normed bounds for 
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the admissible 



< 


k 6 

(iY^ ) can be developed via 

([33[S* (YU)] ) _1 (_S^, ) 

, i i ^i.. ‘ \ 


(4,17) namely 





i i! ([3][B*(Y^)]) _1 j| ; (4.37) 

k s t 

such that c defines the upper bound on j] (^§ m+ i ) !l • Equation 

k 0 

(4.37) defines a pointwise bound on || (iY m+ ^ ) j|. For general use, 
volume averaged values are perhaps more useful. In this context, 

x u 

averaging (4.37) over say the e tn element yields the bounding expression 
!l (i.m + i) e H J V e ]| ( /[D][B*(Y^)] dv)- 1 1] c - (4.38) 


Lastly, we shall develop constraints based on controlling successive 
increments of energy. For the present purposes, this will be considered 
on an element level. In particular, based on (4.8) and Fig. 4.3, the 
energy increment associated with the e in element for the (n+l) tn itera- 
tion is given by the following approximating trapezoid, that is 


F ek 

: m+l 




(F 


m+1 


i e > T 




(4.39) 


such that 
.k 




e „ 


(F k ) e 


[k e ](4Y k tl ) e 


(4.40) 
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In a normed sense, (4.39) can be recast in the following quadratic 
polynomial form in £Y m+ .| !| namely 



k 6 

Solving for |j (aY^ ) jj , we obtain the following element level con- 
straint 



(4.42) 

0 k 

such that here E m+ ^ defines the constraint on the allowable local 
level energy excursion. 

As can be seen from the above development, the bounds on the 
various field variables ultimately resulted in constraints on the 
levels of the aY hierarchy of groups. This is directly due to the 
displacement type of FE formulation used to set up the governing 
equations. 

The foregoing constraints have various advantages depending on 
structural geometry, material type, boundary, and loading conditions. 
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For instance, beam, plate and shell type thin walled structures undergo 
defonnation processes which usually involve large deformation and 
slope changes when over loaded. Such behavior is typically well de- 
fined and thus, can be upper bounded for various loading conditions. 
Once the cotential range of deflection and slope excursion are iden- 
tified, the localized incremental constraints defined on nodal defor- 
mations and displacement gradients can be estimated. In this context, 
and AT . can be obtained directly from design requirements and em- 

** J 

pirical data etc. This is incontrast to stress, strain and energy 
excursions which are less well defined in such structures. Regardless, 
as will be seen later, quite liberal constraints can be employed incon- 
junction with say the AHECS and still maintain stable and efficient 
iteration processes. 

For thick walled structures which are dominated mainly by material 
nonlinearity (plasticity), the use of incremental stress, strain and 
energy excursion controls are more natural. This follows from the fact 
that formally, nonlinear constitutive relations are usually dependent 
on stress, strain and energy excursions. 

Regardless of the approach taken, due to the use of the displace- 
ment type FE formulation, all the foregoing bounds eventually evolve 
into a displacement type constraints In this context, the main thrust 
of the forthcoming benchmarking will be to document the numerical ope- 
rating characteristics of displacement type localized bounds applied 
to the globally constrained INR. Here, special emphasis will be given 
to considering the AHECS type of scheme. 

V. NUMERICAL BENCHMARKING 

In the preceding sections, the constrained INR strategy was 
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modified to include several major improvements. These consist 
of: 

i) The use of localized constraints to help automate/control 
and numerically stabilize the scheme and; 

ii) The use of automatic abscissa/ordinete updating to enable 
the handling of problems which require either load, de- 
flection or mixed load/deflection control. 

To illustrate the improved convergence and efficiency characteristics 
of the modified constrained I NR , the main thrust of this section will 
be to consider the results of several benchmark examples. Such bench- 
marks will attempt to: 

i) Verify the schemes ability to handle highly nonlinear 
problems involving softening/hardening and buckling 
behavior along with inelastic/nonl inear material pro- 
perties and; 

ii) Carefully compare the scheme with the numerical properties 
of the MINR, FINR, and modi fied/BFGS updated INR and HECS 
algorithms. 

Based on such goals, three benchmark examples are "included for 
demonstration purposes namely: 

i) The elastic buckling of a shallow arch; 

ii) The elastic-plastic buckling of shallow arch; 

iii) The load deflection characteristics of a shallow 
spherical cap. 

The choice of the elastic shallow arcli problem depicted 

in Fig. 5-1 follows from its turning point behavior which marks the 

transition between positive and indefinite stiffness behavior. The 
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elastic-plastic behavior of the arch was employed since plasticity 
greatly accentuates the nonlinearity/transitional behavior in the 
vicinity of buckling. Lastly, the spherical cap shown in Fig. 5-2 
was considered since it possesses very shallow stiffness characteris- 
tics followed by strong stiffening. Such properties typically cause 
significant difficulties for unconstrained INR type schemes. 

To implement such goals, the modified constrained INR strategy 
was introduced into the ADINA 77 FE code. This facilitated treat- 
ment of a wide variety of boundary conditions, material types, geo- 
metries and element families. 

As was noted earlier, load stepping in a constrained scheme 
does not converge directly to the designated target values. To 
enable the handling of problems involving specific target loads/ 
deflections, automatic ordinate/abscissa controls must be implemented 
in the constrained methodology. Table 5-1 illustrates a comparison 
of the targeting capabilities of the AHECS and HECS controlled INR. 

As can be seen, the AHECS converges to the specified load increment. 

This is in contrast to the straight HECS whose converged value depends 
on the intersection between the load deflection response and the type 
of constraint surface employed. 

As the first benchmarking activity, we shall consider the capa- 
bility of the so-called AHECS modified INR to handle specific target 
loading values. Of particular interest will be the evaluation of the 
new schemes stability inconjunction with its numerical efficiency. 

For example. Tables 5.2 - 5.4 illustrate the load step sensitivity 
of the AHECS modified INR. As can be seen, with increasing load step 
size, the unconstrained scheme completely fails while the AHECS especially 
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with BFGS updating remains quite stable. Comparing the HECS with the 
AHECS, we see from Tables 5.3 and 5.4 that the number of required 
stiffness updates is significantly reduced with the use of automatic 
abscissa/ordinate adjustment. 

Next, the ability to handle lead stepping in the neighborhood of 
buckling will be considered. Noting the arch depicted in Fig. 5-1, 

Fig. 5-3 illustrates the load deflection characteristics of the % cen- 
trally loaded case. In the context of Table 5.5, the inherent stability 
and efficiency of the AHECS to obtain the turning point is clearly 
illustrated. Again the 5FGS updated version is numerically more effi- 
cient. This is especially true as load increment size increases. • 

Comparison of the handling characteristics in the elastic-plastic 
postbuckling range of arch behavior is given in Table 5.6. Based on 
the model employed, the progressive growth of plastic sites during 
collapse can be seen by noting the sequence defined by Figs. 5.4 and 
5.5. Note the case considered consists of an arch wherein plastici- 
zation is initiated just prior to buckling. This is clearly illustrated 
in Fig. 5.4. Such a zone represents a relatively severe transition in 
behavior and hence a natural benchmark test. As seen from Table 5.6, a 
comparison of the various constraint schemes illustrates that the AHECS 
requires significantly less stiffness updates to yield numerically stable 
converged solution. 

Noting the spherical cap given in Fig. 5.2, Fig. 5.6 illustrates its 
force deflection response to central loading. In the context of such 
system behavior. Table 5.7 depicts the efficiency of the AHECS over the 
straight HECS. Specifically, the AHECS required significantly less full 
stiffness updates to converge to the designated load level. This is 
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particularly true in the shallow range of cap behavior where numerous 
load steps are usually necessary. 

We shall now consider the relative merits associated with the use 
of local constraints to automate load stepping. This will be approached 
from two points of view namely: 

i) The use of local constraints to control stepsize so as to 
yield prescribed loads/displacements and; 

ii) Defining the load/displacement increment sensitivity of 
the scheme. 

For example in terms of the centrally loaded arch. Tables 5.8 - 
5.11 illustrate the aspects denoted by i) and ii) given above. In 
particular, as can be seen from these tables, a wide range of load steps 
and local deflection constraints are stabily and efficiently accomodated 
by the scheme. Specifically, essentially unlimited load step sizes 
can be handled. In addition to illustrating sensitivity to load and 
local constraint step size, the results given in Tables 5.8 - 5.11 also 
point to the fact that local ized constraints can be used to handle either 
load or deflection control. In particular, two target deflection states 
were considered in the postbuckling range namely .4 and .8 inches. These 
are illustrated in Fig. 5-7. As seen from the Tables, the results were 
essentially independent of step size. This is clearly illustrated by the 
fact that the number of iterations does not vary significantly, as step 
size increased. Again, the BFGS version has a decided advantage over the 
MINR type constrained scheme especially as the number of increments is 
increased. 

The preceding developments and associated benchmarks have illu- 
strated the wide ranging potentials/capabilities and numercal robustness 
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of the concept of using localized constraints to, resize and reshape 
the global constraint surface controlling successive iterations of 
the INR family of nonlinear FE solution algorithms. Due to their 
overall algorithmic structure, such procedures are easily adapted to 

use in general puroose FE computer codes. 

VI. SUMMARY 

As seen in the preceding section, this paper has developed 
numerical strategies which enable the automatic adjustment of load/ 
deflection incrementation for NR type nonlinear FE equation solvers. 

This was achieved through the introduction of a hierarchy of localized 
constraints which enable the control/delimiting of the iterative 
excursions of the governing field variables. Such constraints enable 
the automatic adjustment of the size, shape, and orientation of the glo- 
bal constraint surface used to control the NR solver. As can be seen 
from the development, localized bounds can be established for any of 
the governing field variables. This includes the establishment of 
localized bounds on excursions in displacement, the deformation gra- 
dients stress, strain, as well as in energy. 

Due to the generality of the approach taken, a wide variety of 
constraint surface geometries can be handled by the overall methodology. 
Note, such constraints have advantages which are dependent on structural 
configuration, material properties, as well as, on boundary and loading 
conditions. Specifically, deformation and deformation gradient control 
are better suited to thin walled structure. In contrast, for thick 
walled configurations, constraints on stress, strain and energy are more 
useful. For structure containing both thick and thin walled components, 


30 



a combination of localized constraints can be used to control the resizing, 
shaping and orientation process leading to automatic incrementation. 

Due to the manner of formulation, the overall scheme can be 
incorporated in most general purpose FE codes. Such a process would 
require little rearchitecturing, especially if good data management 
procedures are in effect. 

Overall, the project has resulted in numerous papers in a wide 
variety of national and international journals. These include: 

1. Computers and Structures 

2. Computational Mechanics 

3. International Jr. of Engineering Science 

4. Franklin Institute 

5. Finite Elements in Analysis and Design 

6. Engineering Mech. ASCLE 

7. AIAA 

Atotal list of papers is given in Appendix 3. 
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FIGURE AND TABLE CAPTIONS 
Caption 

Intersection of INR type stiffness extrapolation 
of solution curve and F = constant constraint 
curve in force deflection space 

INR type iterations/extrapolations bound by % 
oblique linear, ci rcular' ell ipti c and piecewise 
continuous hyperbolic constraint curves 

Elliptically constrained INR iterations about 
buckl ing point 

Typical failures of oblique linear, circular, 
elliptic and piecewise hyperbolic constraint surfaces 
due to "mis-sizing" 

Load step utilization of various constrained 
INR methodologies 

Rotated and rotated-translated hyperplanar 
constraint 

Iteration process associated with translated 
hyperbolic constraint 

Iteration process associated with abscissa 
stretched HECS 

Targeting oblique 1 inear constraint curve via 
sweeping through fixed angular increments 

Targeting oblique linear constraint curve via 
sweeping through fixed force increments 

Targeting oblique linear constraint curve via 
sweeping through fixed MINR estimated aY 
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Figure No. 


CaDtion 


3.7 


3.3 


4.1 

4.2 

4.3 

5.1 

5.2 

5.3 


5.4 


5.5 


5.6 


5.7 


Geometry associated with successive iterations 
of fixed MI N R estimated *Y sweep of oblique 
1 inear constraint 

Geometry associated with successive iterations 
of fixed MI NR estimated aY stretched elliptic 
constraint 

Geometry of HECS 

Ordinate expansion of HECS 

Incremental energy stored in e tn element 

Arch material properties, geometry and 
boundary conditions 

Spherical cap material properties geometry 
and boundary conditions 

Force deflection behavior of centrally loaded 
arch (Fig. 5.1) 

Elastic-plastic postbuckl ing behavior of 
centrally loaded arch 

Spread of plasticity in postbuckl ing range of 
centrally loaded arch response (Fig. 5.4) 

Load deflection response of centrally loaded 
spherical cap 

Prescribed deflections for centrally loaded 
arch probl em 
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Comparison of single step converged final 
solution values for centrally loaded arch 

Load step sensitivities of FINR, MINR, MINR-AHECS 
and BEGS- AHECS 

Stiffness update/iteration requirements for MINR 
modified HECS and AHECS algorithms: arch subject 

to 30 lb. target load 

Stiffness update/iteration requirements for 
BEGS modified HECS and AHECS algorithms: arch 

subject to 30 lb. target load 

Stiffness update/iteration requirements to 
calculate buckling point of centrally loaded 
arch 

Stiffness update/iteration requirements to 
calculate elastic-plastic postbuckling behavior 
of arch (Fig. 5.4) 

Stiffness update/iteration requirements for cap 
subject to 30 lb. target load 

Comparison of efficiencies of locally constrained 
MINR-HECS and BFGS-HECS for arch with .4 inch 
designated target deflection and 100 lb. aF 

Comparison of efficiencies of locally constrained 
MINR-HECS and BFGS-HECS for arch with .4 inch de- 
signated target deflection and 500 lb. aF 

Comparison of efficiencies of locally constrained 
MINR-HECS and BFGS-HECS for arch with .8 inch 
designated target deflection and 100 lb. a F 

Comparison of efficiencies of locally constrained 
MINR-HECS and BFGS-HECS for arch with .8 inch 
designated target deflection and 500 lb. aF 



Applied Load Step 

, r 

:hecs 

EFGS) 

AHECS 

(*<INR, BEGS ) 

15 lb 

10.35 lb 

15 lb 

25 lb 

16.25 *lb 

25 lb 

30 lb 

19.2 lb 

30 lb 


Table 5.1 Comparison Of Single Step Converged Final Solution 
Values for Centrally Loaded Arch 
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AHECS AHECS 

Applied Load FI NR MI NR (MINR) (BFGS) 

Step Iter. Iter. Iter. Iter 

lb. 



SO Failed Failed 17 6 


Table 5.2 Load Step Sensitivities of FINR, MINR, MINR-AHECS 
and BFGS-AHECS 
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Table 5.3 Stiffness Update/ 1 teration Requirements for Mt NR Modified 
DECS and AHECS Alqori thins: Arch Subject to 30 lb Target. 

































LOAD STEPS 
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Table 5.4 Stiffness Update/ 1 teration Requirements for REGS Modified 
HECS and AIIECS Algorithms: Arch Subject to 30 lb Target 

Load 
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Table 5.5 Stiffness Update/Iteration Requirements to Calculate 
Buckling Point of Centrally Loaded Arch 












































Applied Load Steo 



(lb) 


MINR-HECS 


75 
9 T 


BFGS-HECS 


75 

91 


68 

79 


48 

54 


68 

78 


48 

53 


MINR-AHECS 


51 

no 


4 1 

91 


30 

70 


BFGS-AHECS 


51 

110 


41 

83 


30 

61 


Table 5.6 Stiffness Update/Iteration Requirements to 

Calculate Elastic-Plastic Postbuckl ing Behavior 
of Arch (Fig. 5.4) 
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Table 5.7 Stiffness Update/Iteration Requirements for 
Cap Subject to 30 lb Target Load 









Table 5.8 Comparison of Efficiencies of Locally Constrained 
MINR-HECS and BFGS-HECS for Arch With .4 Inch 
Designated Target Deflection and 100 lb aF 
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No. of Load Steps 

5 

10 • 

25 

50 

100 

Load Step 
(lb) 

500 

500 

500 

500 

500 

6 crit 

(in) 

.08' 

.04 

.016 

.008 

.004 

Final Load Step 
(lb) 

36.3 

36.3 

36.3 

36.3 

36.3 

! 

Final Deflection 
(in) 

.4 

.4 

.4 

.4 

.4 

MNR-HECS 

(Iter) 

141 

110 

99 

91 

129 

BFGS-HECS 

(Iter) 

78 

96 

86 

79 

117 

Table 5.9 Comparison of Efficiencies of Locally Constrained 
MINR-HECS and BFGS-HECS for Arch With .4 Inch 


Designated Target Deflection and 500 lb iF 
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No. of Load Steps 5 10 25 50 TOO 


Load Step TOO TOO TOO TOO 100 

(lb) ‘ 

5 .. .16 .08 .032 .016 .008 

cri t 

(in) . 


Final Load Step 28.6 28.6 28.6 28.6 28.6 

(lb) 


Final Deflection .8 .8 .8 .8 .8 

(in) 


MINR-HECS 180 152 131 128 148 

(Iter) 


BFGS-HECS 90 123 119 112 119 

(Iter) 


Table 5.10 Comparison of Efficiencies of Locally Constrained 

MINR-HECS and BFGS-HECS for Arch with .8m Designated 
Target Deflection and 100 lb aF 
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No. of Load Steps 

5 

10 ■ 

25 

50 

100 

Load Step 
(lb) 

500 

500 

500 

500 

500 

5 crit 

(in) 

.16 

.08 

.032 

.016 

.008 

Final Load Step 
(lb) 

28.6 

28.6 

28.6 

28.6 

28.6 

Final Deflection 
(in) 

.8 

.8 

.8 

.8 

.8 

MNR-HECS 

(Iter) 

191 

162 

141 

139 

157 

MNR-HECS 

(Iter) 

92 

138 

121 

119 

125 


Table 5.11 Comparison of Efficiencies of Locally Constrained 

MINR-HECS and BFGS-HECS for Arch With .8m Designated 
Target Deflection and 500 lb iF 
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fig 3.6 
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fig 4.1 
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